%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Paper: The Technical Default Spread, by Emilio Bisetti, Kai Li, and Jun Yu
% This code is used to produce Table 1, and Column 3 and 4 of Table 2
% Auhtor: Emilio Bisetti, Kai Li, and Jun Yu
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% For dynare package, you need to go to https://www.dynare.org/, download
% the package and install it on your computer. We are using Dynare 4.5.3
% for our model computation.

% If you can not directly run the whole file, you can try to run it section
% by section
%% Produce Table 1 of the papger
clear;
clc;

% Choose the case with risk shocks
cd([cd '/Model with Risk Shocks']);

% Choose the model with covenant
cd([cd '/GE Model computation with covenant']);

load para_v1.mat;

% Following parameters are listed in Table 1
apara_bench=[alppha,gamma,ppsi,beta,delta,zeta,chi,eta,sigma_0,sig1,Ass,xi,rho_a,sigma_a,rho_s0,sigma_s0,sigma_s1]'


%% Load the calibrated parameters
% Aggregation and summary statistics
% For the model with covenant
clear;
clc;

% run dynare code first
dynare mainmodel_tfpunc
% For dynare package, you need to go to https://www.dynare.org/, download
% the package and install it on your computer. We are using Dynare 4.5.3
% for our model computation. 


load para_v1.mat
%load('mainmodel_tfpunc_results.mat');
load('mainmodel_tfpunc_results.mat');
load('simulated_series_raw.mat');
%load('mainmodel_finshock_notd_results.mat');

% Show calibration frequency used;
if freq == 4
    disp('Quarterly Calibration')
elseif freq == 12
    disp('Monthly Calibration')
elseif freq == 1
    disp('Annualy Calibration')
end

aaa = simulated_series_raw(:,:,1);


simulate_years = 400; % take 400 years simulation resutls
num_endo = length(M_.endo_names);
endo_list = cell(num_endo, 1);
discard_period = size(oo_.endo_simul, 2)- simulate_years*freq;

Amoments = NaN(15,options_.simul_replic);

for jj = 1:options_.simul_replic
    simulated_data = simulated_series_raw(:, discard_period+1:end,jj);

    for idx = 1:length(M_.endo_names)
        endo_list{idx} = M_.endo_names(idx,:); % define a cell for each variable
    end

    for idx = 1:length(endo_list)
    %     disp(i) 
        eval([endo_list{idx}, '= transpose(simulated_data(idx,:));']);
    end

    % specify a starting value for F (in level) 
    % recover the K with trend
    length_Q = length(dk); % Length of quarters
    K_Q = zeros(length_Q,1); 
    K_Q(1) = 1; % assign initial value to be 1
    for index = 2:length(dk)
       K_Q(index) = K_Q(index-1)*exp(dk(index));  % recover the time series of K_t
    end

    dK_Q = diff(K_Q)./K_Q(1:end-1);
    Y_Q = exp(y).*K_Q;
    C_Q = exp(c).*K_Q;
    I_Q = exp(ii).*K_Q;
    dY_Q = diff(Y_Q)./Y_Q(1:end-1);
    dC_Q = diff(C_Q)./C_Q(1:end-1);
    act_def_Q = act_def; % actual default probability in level
    TDprob_Q =exp(TDprob);%TDprob;

    N = length(dY_Q);

    % Be careful, we need to treat flow variable and stock variable, and or instant variables differently
    % Turn a column into year X freq matrix 
    % For flow varaible, we need to take the summation of quarterly variables
    % to get annual varibles
    % For output Y
    Y_temp = (reshape(Y_Q, freq, []))'; % 
    Y_a = sum(Y_temp,2);                 % for each year, take summation across freq
    
    % For consumption
    C_temp = (reshape(C_Q, freq, []))';
    C_a = sum(C_temp,2);

    % For investment
    I_temp = (reshape(I_Q, freq, []))';
    I_a = sum(I_temp,2);
 
    dY_a = log(Y_a(2:end)./Y_a(1:end-1));
    dC_a = log(C_a(2:end)./C_a(1:end-1));
    dI_a = log(I_a(2:end)./I_a(1:end-1));

    act_def_a = sum((reshape(act_def_Q, freq, [])), 1)'; % actual default probability in level
    %TDprob_a = sum((reshape(TDprob_Q, freq, [])), 1)'; % technical default probability
    TDprob_a = sum((reshape(TDprob_Q, freq, [])), 1)'; % technical default probability
    
    % Annualized investment to output ratio
    EIY = mean(I_a./Y_a);

    % Auto correlation coefficients of annualized growth rate
    [ac1, junk] = autocorr(dY_a, 1);
    ar1_dy = ac1(2);
    
    % Auto correlation coefficients of annualized consumption growth rate
    [ac1, junk] = autocorr(dC_a, 1);
    ar1_dc = ac1(2);
    
    % Auto correlation coefficients of annualized investment growth rate
    [ac1, junk] = autocorr(dI_a, 1);
    ar1_di = ac1(2);
    
    % leverage regression
    rd_Q = exp(rd);
    rl_Q = exp(rl);
    
    % return: sum up quarterly return to annual ones
    rd_a = sum((reshape(rd, freq, [])), 1)';
    rl_a = sum((reshape(rl, freq, [])), 1)';
    re_a = sum((reshape(re, freq, [])), 1)';
   
    % corelation between consumption growth and investment growth
    COVCI =corr(dI_a, dC_a);
    
    % corelation between consumption growth and investment growth
    COVCY =corr(dY_a, dC_a);
    
    % TD prob and actual default prob
    TDprob_mean = mean(TDprob_a);
    Actual_def_mean = mean(act_def_a);

    amoments_list_PAPER = [100*std(dY_a), 100*std(dC_a), 100*std(dI_a),COVCI,COVCY,ar1_dc, ar1_di,ar1_dy,EIY,mean(exp(rd_a))*100-100,100*std(exp(rd_a)),mean(exp(re_a)-exp(rd_a))*100,mean(exp(rl_a)-exp(rd_a))*100,TDprob_mean*100,Actual_def_mean*100]';

    Amoments(:,jj)=amoments_list_PAPER;
    
end


% produce column 3 of Table 2
Paper_moments = mean(Amoments,2);


fprintf('\n######### Column 4 of Table 2 in the Paper ######### \n');

fprintf('Std of output growth: %.2f \n', Paper_moments(1));
fprintf('Std of consumption growth: %.2f \n', Paper_moments(2));
fprintf('Std of investment growth: %.2f \n', Paper_moments(3));
fprintf('Correlation between consumption and investment growth: %.2f \n', Paper_moments(4));
fprintf('Correlation between consumption and output growth: %.2f \n', Paper_moments(5));
fprintf('Auto-Correlation of consumption growth: %.2f \n', Paper_moments(6));
fprintf('Auto-Correlation of investment growth: %.2f \n', Paper_moments(7));
fprintf('Auto-Correlation of output growth: %.2f \n', Paper_moments(8));
fprintf('Investment to output ratio: %.2f \n', Paper_moments(9));
fprintf('Risk-free rate: %.2f \n', Paper_moments(10));
fprintf('Volaitlity of risk-free rate: %.2f \n', Paper_moments(11));
fprintf('Market return: %.2f \n', Paper_moments(12));
fprintf('Loan spread: %.2f \n', Paper_moments(13));
fprintf('Technical default probability: %.2f \n', Paper_moments(14));
fprintf('Actual default probability: %.2f \n', Paper_moments(15));

%% For the model without covenant
%//////////////////////////////////////////////////////////////////////////

clear;
clc;

cd ..

cd([cd '/GE Model computation nocovenant']);

dynare mainmodel_tfpunc_nocov
% You need to run the dynare code first
% For dynare package, you need to go to https://www.dynare.org/, download
% the package and install it on your computer. We are using Dynare 4.5.3
% for our model computation. 

clear;
clc;

load para_v1.mat % this is the parameter used in the model computation

load('mainmodel_tfpunc_nocov_results.mat'); % this is results generate from runing the dynare code
load('simulated_series_raw.mat'); 

freq = 4;
% Show calibration frequency used;
if freq == 4
    disp('Quarterly Calibration')
elseif freq == 12
    disp('Monthly Calibration')
elseif freq == 1
    disp('Annualy Calibration')
end

simulate_years = 400; % take 400 years simulation resutls
num_endo = length(M_.endo_names);
endo_list = cell(num_endo, 1);
discard_period = size(oo_.endo_simul, 2)- simulate_years*freq;

simulated_data = oo_.endo_simul(:, discard_period+1:end);
Amoments = NaN(14,options_.simul_replic);

for jj = 1:options_.simul_replic
    simulated_data = simulated_series_raw(:, discard_period+1:end,jj);

    for idx = 1:length(M_.endo_names)
        endo_list{idx} = M_.endo_names(idx,:); % define a cell for each variable
    end

    for idx = 1:length(endo_list)
    %     disp(i) 
        eval([endo_list{idx}, '= transpose(simulated_data(idx,:));']);
    end

    % specify a starting value for F (in level) 
    % recover the K with trend
    length_Q = length(dk); % Length of quarters
    K_Q = zeros(length_Q,1); 
    K_Q(1) = 1; % assign initial value to be 1
    for index = 2:length(dk)
       K_Q(index) = K_Q(index-1)*exp(dk(index));  % recover the time series of K_t
    end

    dK_Q = diff(K_Q)./K_Q(1:end-1);
    Y_Q = exp(y).*K_Q;
    C_Q = exp(c).*K_Q;
    I_Q = exp(ii).*K_Q;
    act_def_Q = act_def; % actual default probability in level

    % Be careful, we need to treat flow variable and stock variable, and or instant variables differently
    % Turn a column into year X freq matrix 
    % For flow varaible, we need to take the summation of quarterly variables
    % to get annual varibles
    % For output Y
    Y_temp = (reshape(Y_Q, freq, []))'; % 
    Y_a = sum(Y_temp,2);                 % for each year, take summation across freq
    % For consumption
    C_temp = (reshape(C_Q, freq, []))';
    C_a = sum(C_temp,2);
    % For Investment
    I_temp = (reshape(I_Q, freq, []))';
    I_a = sum(I_temp,2);

    dY_a = log(Y_a(2:end)./Y_a(1:end-1));
    dC_a = log(C_a(2:end)./C_a(1:end-1));
    dI_a = log(I_a(2:end)./I_a(1:end-1));
    act_def_a = sum((reshape(act_def_Q, freq, [])), 1)'; % actual default probability in level

    % Annualized Investment to output ratio
    EIY = mean(I_a./Y_a);

    % Auto correlation coefficients of annualized growth rate
    [ac1, junk] = autocorr(dY_a, 1);
    ar1_dy = ac1(2);
    
    % Auto correlation coefficients of annualized consumption growth rate
    [ac1, junk] = autocorr(dC_a, 1);
    ar1_dc = ac1(2);
    
    % Auto correlation coefficients of annualized investment growth rate
    [ac1, junk] = autocorr(dI_a, 1);
    ar1_di = ac1(2);
    
    % return: sum up quarterly return to annual ones
    rd_a = sum((reshape(rd, freq, [])), 1)';
    rl_a = sum((reshape(rl, freq, [])), 1)';
    re_a = sum((reshape(re, freq, [])), 1)';

    % corelation between consumption growth and investment growth
    COVCI =corr(dI_a, dC_a);

    % corelation between consumption growth and investment growth
    COVCY =corr(dY_a, dC_a);

    %TDprob_mean = mean(TDprob_a);
    Actual_def_mean = mean(act_def_a);
    amoments_list_PAPER = [100*std(dY_a), 100*std(dC_a), 100*std(dI_a),COVCI,COVCY,ar1_dc, ar1_di,ar1_dy,EIY,mean(exp(rd_a))*100-100,100*std(exp(rd_a)),mean(exp(re_a)-exp(rd_a))*100,mean(exp(rl_a)-exp(rd_a))*100,Actual_def_mean*100]';
    Amoments(:,jj)=amoments_list_PAPER;
end

% calucatethe mean of all replications
Paper_moments = mean(Amoments,2);

fprintf('\n######### Column 4 of Table 2 in the Paper ######### \n');

fprintf('Std of output growth: %.2f \n', Paper_moments(1));
fprintf('Std of consumption growth: %.2f \n', Paper_moments(2));
fprintf('Std of investment growth: %.2f \n', Paper_moments(3));
fprintf('Correlation between consumption and investment growth: %.2f \n', Paper_moments(4));
fprintf('Correlation between consumption and output growth: %.2f \n', Paper_moments(5));
fprintf('Auto-Correlation of consumption growth: %.2f \n', Paper_moments(6));
fprintf('Auto-Correlation of investment growth: %.2f \n', Paper_moments(7));
fprintf('Auto-Correlation of output growth: %.2f \n', Paper_moments(8));
fprintf('Investment to output ratio: %.2f \n', Paper_moments(9));
fprintf('Risk-free rate: %.2f \n', Paper_moments(10));
fprintf('Volaitlity of risk-free rate: %.2f \n', Paper_moments(11));
fprintf('Market return: %.2f \n', Paper_moments(12));
fprintf('Loan spread: %.2f \n', Paper_moments(13));
fprintf('Actual default probability: %.2f \n', Paper_moments(14));
